{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import string\n",
    "#https://stackoverflow.com/questions/19726663/how-to-save-the-pandas-dataframe-series-data-as-a-figure\n",
    "import six"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "tbpath = \"../../fits/\"\n",
    "productpath = \"../../postfit_derivatives/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "models = [\"fulllinearmodel_fit_table.csv\",\"reducedlinearmodelNegBinom_fit_table.csv\",\n",
    "          \"reducedlinearmodelq0_fit_table.csv\",\"reducedlinearmodelq0ctime_fit_table.csv\",\n",
    "         \"nonlinearmodelq0ctime_fit_table.csv\",\"nonlinearmodel_fit_table.csv\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "rois = []\n",
    "for model in models:\n",
    "    df = pd.read_csv(tbpath + model) #get rois in all tables (some may have failed)\n",
    "    rois += list(df.roi.unique())\n",
    "    \n",
    "rois = list(set(rois))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gambia\n"
     ]
    }
   ],
   "source": [
    "dfreport = pd.DataFrame(columns=['Region', 'Model','median IFR0','CI lower','CI upper','median IFRc','CI lower','CI upper','delta weeks'])\n",
    "k = -1\n",
    "for roi in rois:\n",
    "    k += 1\n",
    "    model = \"fulllinearmodel_fit_table.csv\" \n",
    "    df = pd.read_csv(tbpath + model)\n",
    "    theta = \"ifr (week 0)\"\n",
    "    try:\n",
    "        measure1 = df.loc[(df.roi==roi)&(df['quantile']==0.5),theta].values[0]\n",
    "        measure2 = df.loc[(df.roi==roi)&(df['quantile']==0.025),theta].values[0]\n",
    "        measure3 = df.loc[(df.roi==roi)&(df['quantile']==0.975),theta].values[0]\n",
    "        x = [roi,model.split('_fit_table.csv')[0],np.round(measure1,4),np.round(measure2,4),np.round(measure3,4)]\n",
    "        #find latest week with data\n",
    "    #     maxweek = 11\n",
    "        for i in np.arange(11,0,-1):\n",
    "            theta = 'ifr (week '+str(i)+')'\n",
    "            x2 = df.loc[(df.roi==roi)&(df['quantile']==0.5),theta].values[0]\n",
    "            if np.isfinite(x2):\n",
    "    #             print(x2)\n",
    "                measure1 = df.loc[(df.roi==roi)&(df['quantile']==0.5),theta].values[0]\n",
    "                measure2 = df.loc[(df.roi==roi)&(df['quantile']==0.025),theta].values[0]\n",
    "                measure3 = df.loc[(df.roi==roi)&(df['quantile']==0.975),theta].values[0]\n",
    "                x += [np.round(measure1,4),np.round(measure2,4),np.round(measure3,4),i]\n",
    "                break\n",
    "    except:\n",
    "        print(roi)\n",
    "    try:\n",
    "        dfreport.loc[k] = x\n",
    "    except:\n",
    "        print(roi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0049\n",
      "0.025    0.00132\n",
      "0.500    0.00490\n",
      "0.975    0.03832\n",
      "Name: median IFR0, dtype: float64\n",
      "0.0228\n",
      "0.025    0.00190\n",
      "0.500    0.02280\n",
      "0.975    0.13738\n",
      "Name: median IFRc, dtype: float64\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1MAAAE0CAYAAADJ1mlXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd1RU59YG8OfAIB1Fmg1E0YASBUUBG+BYQGNBYq9YiBolRsSCMYIFJZZosF4rWC9eRL1GQAMD3sQumuunYkkATdSgGEEQpM73h4u5nMxQHAtJeH5rzUrY5z3v2WcgZDZvOYJcLgcRERERERG9Ho3aToCIiIiIiOiviMUUERERERGRGlhMERERERERqYHFFBERERERkRpYTBEREREREamBxRQREREREZEaJNUcf7N90z08Xv0zOfmNuqkTIiJe/dPXtzazICIiIiIiMaGyAxyZIiIiIiIiUgOLKSIiIiIiIjWwmCIiIiIiIlIDiykiIiIiIiI1sJgiIiIiIiJSA4spIiIiIiIiNbCYIiIiIiIiUgOLKSIiIiIiIjVU99BeIiIi+hvLyclBVlYWioqKajsVIqL3pl69ejA1NUX9+vXfqB8WU0RERHXUy5cvkZmZiWbNmkFXVxeCINR2SkRE75xcLkdBQQF+/fVXaGtrQ0dHR+2+OM2PiIiojnry5AnMzMygp6fHQoqI6gxBEKCnpwdTU1M8efLkjfp6o5GptIg9VR5v/FsmAED3TS5CRERE78TLly/RqFGj2k6DiKhWGBoa4unTp2/UB0emiIiI6qiSkhJIJJzxT0R1k0QiQUlJyRv1wWKKiIioDuP0PiKqq97G7z8WU0RERERERGpgMUVERET0niQnJ0MQBERERFQZo1cEQYCvr29tp0FUKRZTRERERPS3dfv2bXh7e8PY2Bj6+vro0aMHZDJZbadFfxNcdUpERERUi9zc3FBQUAAtLa3aTuVv5+eff0bXrl0hkUgwb9481K9fH9u3b4enpyfi4uLQu3fv2k6R/uJYTBERERHVIg0NjTd6aChVLigoCNnZ2UhJSYGjoyMAYPz48bC3t8eMGTNw69YtbsJCb4TT/IiIiOhvLSIiAoIgIDExEUuXLkXz5s2hq6sLFxcXnD9/HgBw+vRpdO/eHfr6+mjcuDGWLVumsq/Lly9jyJAhMDU1hba2NmxtbREaGqpye+Vjx46hQ4cO0NHRgaWlJb788ksUFxcrtVO1ZqqsrAyhoaFwc3NDo0aNUK9ePVhZWWH69OlKz8XJyMiAIAgICQnBt99+i86dO0NHRweNGzfG3Llza7T1c8+ePWFtbS2KHTx4EIIgwMHBQRTfsmULBEHAhQsXFDG5XI4tW7bAyckJenp6MDAwQM+ePZGUlKTyelFRUejevTsMDQ2hp6cHFxcXREdHV5snAFy5cgWNGjVC27Ztcf/+/UrbvXjxAv/+97/h4eGhKKQAwMDAAFOmTMGdO3dw6dKlGl2TqDIspoiIiKhOWLBgAY4ePYpZs2YhODgYaWlp6Nu3L44ePQofHx/06NEDa9asgZ2dHRYvXox9+/aJzj9x4gS6deuGO3fuYM6cOQgPD0eXLl2wePFijBo1StT2yJEjGDJkCHJycrB48WL4+/sjOjoa8+fPr1GuRUVFWL16NVq3bo25c+ciPDwcffr0wc6dO+Hh4YGioiKlc2JjYzFp0iT069cP69atg4ODA9asWYNVq1ZVez2pVIp79+7h559/VsQSExOhoaGB//u//0NWVpYiLpPJYGRkhE6dOili48aNw8yZM9GqVSusWrUKS5YsQU5ODvr06YN///vfomstWrQII0eOhKGhIZYtW4awsDDo6elh2LBh2LRpU5V5njx5Eu7u7rCxscEPP/wAKyurStteu3YNhYWF6NKli9IxV1dXAGAxRW9OLpdX9arSz7sjq3zl29rJ823tquuG5HK5fPfuVy8iIqL35ObNm7Wdwnuxe/duOQB5hw4d5IWFhYr4sWPH5ADkEolEfunSJUW8sLBQ3qhRI7mrq6siVlBQILewsJD36NFDXlxcLOr/66+/lgOQJyUlyeVyubykpERuaWkpNzExkT958kTRLjs7W25lZSUHIN9d4f/5SUlJSrGysjJ5fn6+0r3s2LFDDkAeFRWliKWnp8sByPX09OTp6emiPuzt7eWNGjWq9j364Ycf5ADk27ZtU8RatGghHzt2rOh6ZWVlclNTU/mAAQMU7WJiYuQA5P/4xz9EfRYXF8udnJzk1tbW8rKyMrlcLpenpKTIAciDgoKUchg8eLDc0NBQ/vz5c0UMgHzChAlyuVwu37Nnj1xLS0s+ePBgle/NH0VHR8sByDdv3qx07MaNG5XmQXVLDX8PVlovcc0UERERiX3+OfDjj7WdhZijI7B+/Rt1MX36dNSrV0/xdY8ePQAALi4uolGWevXqwdnZGWfOnFHEvvvuO2RmZmLlypXIzs4W9du/f38EBATg1KlT8PDwQEpKCn755RcEBgbC1NRU0a5+/fqYNm0aFi5cWG2ugiBAV1cXAFBaWorc3FyUlJRAKpUCAC5cuIDhw4eLzvH29hZN1RMEAT179sTGjRuRl5cHAwODSq/n7OwMAwMDyGQy+Pn54d69e0hPT8fGjRvx448/IjExEcOHD1eMUpXnAQD79u2DoaEhvL29RSNYADBw4ECEhITg7t27+OCDD7B//34IgoAJEyYotR00aBCOHTuGc+fOoW/fvqJjYWFhWLhwIT755BNs2rQJmpqa1b6H+fn5AABtbW2lY+Vr1MrbEKmLxRQRERHVCS1bthR9bWxsDABo0aKFUltjY2PR2qTU1FQAwKRJkyrtPzMzEwCQlpYGALCzs1Nq07Zt2xrne+jQIaxduxZXr15VWmv17NkzpfZ/vD8AMDExAQA8ffq0ymJKS0sL3bt3V6xxSkxMhEQigZubG6RSKWJjYwFAsaV4xWIqNTUVubm5sLCwqLT/zMxMfPDBB0hNTYVcLlf53lRsW1FMTAxyc3Ph5+eHrVu3VnreH+np6QEACgsLlY69fPlS1IZIXSymiIiISOwNR4D+rCobzajJKIdcLgcArF69WrSZQUVNmjRRP7k/iImJwYgRI+Ds7IxvvvkGlpaW0NHRQWlpKby8vFBWVqZ0TlX3UZ5/VaRSKeLj43Hjxg3IZDJ07twZBgYGkEqlCA8Px/379yGTyWBqaor27duL+jYzM8OBAwcq7fvDDz9UtBUEAXFxcZXma29vL/ra2dkZGRkZiI6OxieffCIaRaxK+ffjwYMHSsfKY02bNq1RX0SVYTFFREREVI3WrVsDAPT19at9NlH5CNGtW7eUjt28ebNG19u7dy90dHSQlJQkGj1R1efbUj7alJiYCJlMhsmTJwMAPDw8oKmpiVOnTuE///kP+vTpI9pOvHXr1rhz5w5cXV2rHP0qbxsfHw8rKyu0adOmRnk1a9YMkZGRkEql6N27N+Lj4xUbSFSlXbt20NbWxrlz55SOle/iWNPCjKgy3M2PiIiIqBqenp4wNzdHWFgYfv/9d6XjBQUFyM3NBQA4OTmhWbNm2L17t2hd0PPnz2s8TU1TUxOCIIhGoORyOZYvX/6Gd1K5Dh06wNjYGFu3bsWjR48UxVX9+vXRsWNHrFu3Djk5OaIpfsCr5zaVlZUhKChIZb8Vp+2NGzcOALBw4UKUlpZW2baipk2b4vTp02jSpAn69u0rWs9WGQMDAwwcOBDJycn473//q4jn5eVhx44daN26NZydnavth6gqHJkiIiIiqoa+vj727NkDb29v2NraYtKkSWjVqhWys7Nx69YtxMTE4MiRI4pRnHXr1mH48OFwdnaGn58fJBIJdu3aBRMTkyqfjVRu6NChOHz4MKRSKcaPH4/i4mIcPXr0nW6YoKGhAXd3dxw9ehQ6Ojro2rWr4phUKsVXX32l+Pc/5jpx4kRs3LgRV65cwYABA2Bqaopff/0V586dw08//aRYR9a5c2eEhIQgJCQEjo6OGDZsGJo0aYJHjx4hJSUFsbGxKrd9B4BGjRohOTkZvXv3hpeXF7799lu4u7tXeU8rV65EYmIi+vbti9mzZ8PIyAjbt2/HgwcPcOLECT6wl94YR6aIiIiIasDT0xOXLl2Cp6cn9u3bhxkzZmDNmjVITU1FQECAaB3R0KFDER0dDSMjI4SEhCA8PBxDhw5VFCTVGTlyJLZt24a8vDwEBgZi1apVsLW1xcmTJ9/V7QH4X6HUtWtX0S54vXr1AvBqhMjW1lbpvF27dmHPnj3Q0NDAypUr4e/vj8jISBgYGGDlypWitsHBwfj222/RpEkTrF+/HjNmzMC2bdtQWFiI8PDwKvMzNzdHUlISWrVqhf79+yMxMbHK9q1atcKZM2fg6uqKsLAwBAYGQl9fH/Hx8fD09KzRe0JUFaGaBYlVHkyL2FNl543DXv3Ho3sr9XXzqnvKn3ru61ubWRARUR2Smppa43UrRER/RzX8PVjpECZHpoiIiIiIiNTAYoqIiIiIiEgNLKaIiIiIiIjUwGKKiIiIiIhIDSymiIiIiIiI1MBiioiIiIiISA0spoiIiIiIiNTAYoqIiIiIiEgNLKaIiIiIiIjUwGKKiIiIiIhIDSymiIiIiIiI1MBiioiIiIiISA0spoiIiIjek+TkZAiCgIiIiCpjdU1ISAgEQUBGRoYiFhERAUEQkJycXGt5VcfDwwPW1ta1nQbVIhZTRERERCQyYsQICIKAXr161XYqNfLo0SN88cUX8PLygpmZGQRBgK+vb6Xtra2tIQiCyldWVpaobUREBNavX/+O70BMLpdj3759GDlyJFq1agU9PT1YWVlh0KBBuHDhgspzysrKsG7dOtjZ2UFHRweWlpaYM2cOXrx4IWr37NkzfPPNN+jbty8sLS2hq6sLW1tbfPLJJ/jll19U9p2TkwN/f380bdoUOjo6sLe3x5YtWyCXy9/6vf/VSGo7ASIiIqK6zM3NDQUFBdDS0qrtVAAAT58+xdGjR2FjY4OkpCRkZGS889GXRYsWYcGCBdDW1lbr/Nu3b2PFihWwtLRE586dERcXV+05dnZ2+OKLL5TihoaGoq8jIiKQkZGBzz//XK3c1FFYWIhx48bB0dERI0eORIsWLfDo0SNs3boVXbp0wZ49ezB27FjRObNnz0Z4eDiGDBmCOXPmIDU1FeHh4bh69SoSEhKgofFqDOXChQuYM2cOevXqhZkzZ8LU1BTXr1/HP/7xDxw6dAhnz55F27ZtFf0WFRWhT58+uHr1Kvz9/dGmTRvExcXh008/RWZmJkJCQt7b+/JnxGKKiIiIqBZpaGhAR0enttNQ2LdvH4qLixEVFYUuXbpg9+7dWLJkyTu5Vm5uLgwNDSGRSCCRqP+x1MnJCY8fP4aZmRmysrJgZmZW7TkWFhZKBcmfhUQiQXJyMtzd3UVxPz8/2NvbY86cORg9erSiQLpx4wY2bNgAHx8fHD58WNG+RYsW+Oyzz/DPf/4To0ePBvCqiLx9+zZsbGxEfX/00Ufo06cPFi9ejOjoaEV8x44duHTpEsLDw+Hv76/I4+OPP8aKFSswceJENG/e/J28D38FnOZHREREf2vla28SExOxdOlSNG/eHLq6unBxccH58+cBAKdPn0b37t2hr6+Pxo0bY9myZSr7unz5MoYMGQJTU1Noa2vD1tYWoaGhKCkpUWp77NgxdOjQQTHl6ssvv0RxcbFSO1VrpsrKyhAaGgo3Nzc0atQI9erVg5WVFaZPn46nT5+Kzs/IyIAgCAgJCcG3336Lzp07Q0dHB40bN8bcuXNV5laVnTt3wsPDA05OThgwYAAiIiJQVlam1M7X1xeCIODJkycYP348TExMoK+vj169euHKlSuV5hgVFQUnJyfo6uoqPpyrWjP1OgwNDWtUQP1RSUkJnj9/Xulxa2trnD59Gvfu3RNNBfzjOq6HDx9i1KhRMDY2hp6eHjw9PXHnzp3XzqecRCJRKqSAVwWgu7s7Hj9+jMePHyviBw8ehFwuVxo98/Pzg56eHvbt2ye6pz8WUgDQu3dvNGzYENevXxfFDxw4AD09Pfj5+Ynin3/+uaLorstYTBEREVGdsGDBAhw9ehSzZs1CcHAw0tLS0LdvXxw9ehQ+Pj7o0aMH1qxZAzs7OyxevFj0ARQATpw4gW7duuHOnTuYM2cOwsPD0aVLFyxevBijRo0StT1y5AiGDBmCnJwcLF68GP7+/oiOjsb8+fNrlGtRURFWr16N1q1bY+7cuQgPD0efPn0UhU5RUZHSObGxsZg0aRL69euHdevWwcHBAWvWrMGqVatq/B5dunQJ//d//4cJEyYAeFUw3b9/HwkJCZWe4+XlhUePHiEkJASff/45Ll++DHd3d6UP5QBw9OhRTJ8+HV5eXggPD0e/fv1qnNvbduHCBejp6aF+/fpo0KABJkyYgIcPH4rarF+/HnZ2djA1NcXevXsVrzZt2ijavHjxAm5ubtDU1MSKFSswc+ZMJCcnY/DgwSgtLX3ref/666+oV68eGjRooIhdunQJGhoacHZ2FrXV0dGBo6MjLl26VG2/OTk5yM3NhYWFhSJWVlaGK1euKP4oUJGzszMEQahR339nVY6npkXseV95EBEREb1TpaWlOH/+POrVqwcAaNu2LQYPHoxhw4bh3Llz6NSpEwBg8uTJaN68OTZt2qSYBvby5UtMnjwZLi4ukMlkiilpU6dOhYODAwICApCcnAwPDw+UlpZi1qxZaNiwIS5evAhTU1NF2/bt29coV21tbTx69Ai6urqK2LRp09C1a1dMmTIFR48exfDhw0Xn3LhxAzdu3FCsb5o2bRratWuHDRs2YOHChTW67q5du6Cvr4+PP/4YANCvXz+YmZlh586d6Nu3r8pzmjdvjsOHD0MQBACAj48POnfujMDAQMTHxyvleO3aNVExUhvs7e0xZcoUtGnTBsXFxUhOTsaOHTuQmJiIixcvokmTJgAAb29vrF+/HgUFBZVOCczKysLcuXMxb948RczMzAzz5s1DQkICPD0931resbGxuHjxIsaNGycqbh4+fKgYLf2jpk2b4uzZsygqKlL87KsSGhqK4uJiRSENvNqsoqCgAE2bNlVqr62tDVNTUzx48OAN7+qvjWumiIiISCw+Hvjtt9rOQqxRI8DL6426mD59uujDZI8ePQAALi4uikIKAOrVqwdnZ2ecOXNGEfvuu++QmZmJlStXIjs7W9Rv//79ERAQgFOnTsHDwwMpKSn45ZdfEBgYqCikAKB+/fqYNm1ajQobQRAUhVRpaSlyc3NRUlICqVQK4NWoyh+LKW9vb9FGEYIgoGfPnti4cSPy8vJgYGBQ5TULCgpw8OBBfPzxx4q2WlpaGDNmDLZs2YLff/8dDRs2VDpv3rx5ikIKeLV+qU+fPkhISFC67kcffVTrhRTwapSxopEjR8LNzQ1jxoxBcHAwtm/fXuO+NDQ08Nlnn4li5d+nu3fvvrVi6u7duxg3bhyaNm2KtWvXio7l5+dXunlHedGVn59faTEVHR2NNWvWwMvLCxMnThT1C6DKvsvb1FWc5kdERER1QsuWLUVfGxsbA3i1SP+PjI2NRWuTUlNTAQCTJk2CmZmZ6GVnZwcAyMzMBACkpaUBgCJeUcVd0qpz6NAhuLi4QFdXF8bGxjAzM1Pcw7Nnz6q9PwAwMTEBAKV1VqpER0cjJycH7u7u+OmnnxQvNzc3FBYWKk17LKeqOGrbti1KS0tx7949UfyDDz6oNo/aMnr0aFhbWysVWtVp0qSJ0hS413nfayI9PR29evWCIAiIi4tTWh+mp6eHwsJClee+fPlS0UaV2NhYjBkzBk5OToiKihIVxuXnVNV3Zf3WFRyZIiIiIrE3HAH6s9LU1HyteEXlz9NZvXo1HB0dVbYpnxr2NsTExGDEiBFwdnbGN998A0tLS+jo6KC0tBReXl4qN4So6j5q8jygnTt3Ang1zVGVXbt2KY3AvK4/+wdva2tr0YhkTbzp+16djIwM9OzZE3l5eUhMTES7du2U2jRp0gQ3b95EYWGh0ijSgwcPYGpqqnJUKj4+Hj4+PrC3t8epU6dgZGQkOm5sbAxdXV2VU/kKCwuRlZWlcqOMuoTFFBEREVE1WrduDQDQ19dH7969q2xbPkJ069YtpWM3b96s0fX27t0LHR0dJCUliQoQVX2+DT///DP+85//YMyYMfD29lY6npiYiK1btyIlJQVOTk6iY6mpqXB1dRXFbt68CU1Nzb/cltk//fSTaAMGAKKRmvctIyMDHh4eyMnJQUJCAjp06KCyXefOnXHq1ClcvHhRMX0VeDVy9OOPP8LNzU3pnPj4eHh7e8POzg4JCQmKkdqKNDQ00LFjR1y9elWpULt48SLkcrloimxdxGl+RERERNXw9PSEubk5wsLC8PvvvysdLygoQG5uLoBXa4aaNWuG3bt3IysrS9Hm+fPn2Lp1a42up6mpCUEQRCNQcrkcy5cvf8M7UW3Xrl2Qy+UICAjA0KFDlV7luxDu2rVL6dxVq1aJRmCuXLmChIQE9OrVq9p1WrVB1fcPADZt2oRff/0VAwcOFMUNDAzw7NmztzLK9Dru3buHnj17Ijs7G6dOnVIqYisaMWIEBEHA+vXrRfHt27cjPz8fY8aMEcVPnTqFIUOGwNbWFomJiSrXwpUbNWoU8vPzsW3bNlF8/fr1kEgkGDFihBp39/fBkSkiIiKiaujr62PPnj3w9vaGra0tJk2ahFatWiE7Oxu3bt1CTEwMjhw5Ag8PD2hqamLdunUYPnw4nJ2d4efnB4lEgl27dsHExAT379+v9npDhw7F4cOHIZVKMX78eBQXF+Po0aPvZLF/aWkpIiIiYG1tjY4dO6psY21tDScnJxw4cABr164VrRG6d+8ePD09MWjQIDx69AgbN26Erq4uVq9e/dZzrUp5oVn+Hl27dk0Rc3NzU4zO7NmzBzt37oSXlxesra1RUlKC5ORkHD16FDY2NkoPKHZ1dcW3336LmTNnomvXrtDU1IRUKoW5uflr52htbY179+5VW5jl5uaiZ8+eyMjIgL+/P27fvo3bt2+L2vTp00cxitauXTvMmDEDGzduhI+PD/r374/U1FSEh4fD3d1d8cBe4NWz0gYPHgy5XI6JEyciLi5O6foVdy708/PD7t27ERAQgIyMDLRp0waxsbE4cuQIFi1aJNr0pC5iMUVERERUA56enrh06RLCwsKwb98+PHnyBMbGxrCxsUFAQIBo2/OhQ4ciOjoaS5cuRUhICMzNzeHr6ws3N7dKtxivaOTIkcjNzcW6desQGBgIY2NjDBw4EGFhYYrNDd6W+Ph4PHz4EAEBAVW2+/jjj7Fw4ULExMSIPpzHx8cjICAAwcHBKCgogKurK1avXl3jbeDfli+//FL09dWrV3H16lUAQHBwsKKY6ty5M2QyGaKiovDkyRPI5XK0aNEC8+fPx4IFC0TPbwKA2bNnIy0tDdHR0di6dSvKysqQlJSkVjGVl5dXo7V1T58+RXp6OgBgw4YNKtskJSWJpiSuX78e1tbW2LZtG06cOAFTU1P4+/tj6dKl0ND432S069evKzalmD17tsq+KxZT9erVQ0JCAhYtWoSDBw/i6dOnsLGxwYYNGzBjxozqb/pvTqiqMk6L2PNG45mNw1YCAHRvpb5JN3VD+VPPfX1rMwsiIqpDUlNT/xTbVNNfk6+vLyIjI9/79Le/qmvXrsHBwQG7du0SbT9OtauGvwcrXTjHNVNERERERO/YyZMn4eDgIHooLv31sZgiIiIiInrH5s6dix9//FE05Y7++vjdJCIiIiIiUgOLKSIiIiJ6bREREVwvRXUeiykiIiIiIiI1sJgiIiIiIiJSA4spIiIiIiIiNbCYIiIiIiIiUgOLKSIiIiIiIjWwmCIiIiIiIlIDiykiIiIiIiI1sJgiIiIiIiJSA4spIiIioncgOTkZgiAgIiKitlNRW0hICARBQEZGRo3aC4IAX1/fd5rTn8Hf4XtLbweLKSIiIiKiWhAbG4uuXbtCX18fDRs2xLBhw5Cenl7badFrkNR2AkRERPTnkxaxp7ZTEGnpO762U3htbm5uKCgogJaWVm2norZFixZhwYIF0NbWru1U/nZiYmIwdOhQODg4YPXq1cjJycH69evRrVs3XL58GU2aNKntFKkGWEwRERFRnVNaWorCwkLo6em9s2toaGhAR0fnnfX/PkgkEkgk/Lj4thUXF8Pf3x+Wlpb4/vvvYWBgAADo168fnJycEBISgm3bttVyllQTnOZHREREf2sREREQBAEJCQlYtmwZbGxsoKOjg0OHDinaXL58GUOGDIGpqSm0tbVha2uL0NBQlJSUKPV3+PBhODg4QEdHB1ZWVliyZAkSEhKU1tBUtq7mxYsXCAoKgo2NDbS1tdGoUSOMHz8e9+7dE7WreP7u3bthb28PbW1tNG/eHKtWrarRvd+7dw+CICA4OFgU9/T0hCAIWLdunSju4uKCNm3aKL6ubM3UjRs34OXlpZieNmbMGDx+/LjSPKKiotC9e3cYGhpCT08PLi4uiI6Orjb/wsJC6OrqYsKECaL41KlTIQgCZs2aJYqPGDECRkZGou9bTk4O5s+fj1atWkFbWxtmZmYYNWoU0tLSVF5vxYoVsLe3h46ODho0aICBAwfi6tWr1eYKAJGRkdDS0sLQoUPx8uXLStudPn0aDx8+xJQpUxSFFAA4OjrCw8MDUVFRKC4urtE1qXbxTw1ERERUJwQGBqK4uBh+fn4wMjKCra0tAODEiRPw8fFBq1atMGfOHDRs2BDnzp3D4sWL8eOPP+Jf//qXoo+oqCiMGjUKNjY2CA4OhkQiQWRkJI4fP16jHIqLi+Hp6YkzZ85g6NChmDNnDu7evYstW7bg1KlTuHz5Mpo1ayY6Z+vWrcjMzMTkyZPRoEED7Nu3D/Pnz0ezZs0wevToKq/XvHlztGzZEjKZDEuWLAEAFBUV4YcffoCGhgZkMhlmz54NAHj+/DlSUlIwderUKvtMT09Hjx49UFhYiJkzZ8LS0hLHjx+Hl5eXyvaLFi1CaGgovLy8sGzZMmhoaODIkSMYNmwYNm7ciBkzZlR6LW1tbXTt2n4lBT4AACAASURBVBVJSUmieGJioiL/cnK5HMnJyejRo4diNC0nJwddu3bF/fv3MWnSJNjb2+PRo0fYvHkzXFxccPnyZTRv3hzAq++Nl5cXzp49i3HjxmHmzJnIycnB9u3b0a1bN/znP/9Bp06dKs11xYoV+OKLLzBjxgyEh4dDQ6PyMYtLly4BALp06aJ0zNXVFTKZDHfu3IG9vX2lfdCfA4spIiIiqhMKCgpw9epV0dS+ly9fYvLkyXBxcYFMJlN8CJ86dSocHBwQEBCA5ORkeHh4oKSkBAEBATAzM8PFixdhbGwMAJg+fTrat29foxwiIiJw5swZzJ07VzS61Lt3bwwYMABBQUHYu3ev6Jz79+8jNTUV9evXBwBMmjQJzZs3x4YNG6otpgBAKpUiMjIS+fn50NPTw/nz55Gfn4+xY8fi2LFjKCkpgUQiwenTp1FaWgqpVFplf1988QWePXsGmUyGnj17AgBmzJgBHx8fpRGcK1euIDQ0FEFBQVixYoUi/tlnn8Hb2xtBQUEYP348DA0Nq8xfJpPh7t27aN26Ne7fv4+ff/4ZY8eOxb59+5CZmQkLCwtcv34djx8/FuW/ePFipKWl4fz583BwcFDEfX190a5dOwQHBytGDjdu3Ijk5GTEx8fD09NT0fbTTz/Fhx9+iMDAQCQnJyvlV1ZWBn9/f2zevBmhoaFYuHBhle8fADx8+BAA0LRpU6Vj5bEHDx6wmPoL4DQ/IiIiqhOmT5+utEbqu+++Q2ZmJiZOnIjs7GxkZWUpXv379wcAnDp1CgCQkpKChw8fwtfXV1FIAYCBgQGmTZtWoxyOHDkCDQ0NBAUFieIfffQRHB0dcezYMZSVlYmOTZw4UVFIAYCenh5cXV1x9+7dGl1TKpWiuLgY33//PQBAJpPB3Nwcs2bNQm5urmKUJCkpCYIgKAokVcrKynD8+HF06tRJ1E4QBMybN0+p/f79+yEIAiZMmCB6b7OysjBo0CDk5ubi3Llz1eZfnnf5PzU1NRVTEMvj5aNX5e3lcjn2798PNzc3NG3aVHRtfX19uLq6Kr63ALBv3z7Y2dnByclJ1LaoqAh9+vTBDz/8gIKCAlFuL1++xNChQ7Ft2zZERETUqJACgPz8fABQubFH+Tq78jb058aRKSIiIqoTPvjgA6VYamoqgFejPZXJzMwEAMWW1eXTAytSFVMlPT0dTZo0ERVj5ezt7fHjjz8iKysL5ubminjLli2V2pqYmODp06eKr3NycpQ+6JuZmUFTU1NUjHh6eipGlDp27AhjY2PIZDJ06dIFMpkMDg4OaNiwYaX5P378GHl5ebCzs1M61rZtW6VYamoq5HK5yvblyt/fynTu3BmGhoaQyWSYOnUqZDIZOnXqBBsbG7Rr1w4ymQyjRo2CTCZDw4YN4ejoCAB48uQJnj59ilOnTsHMzExl3xWn4qWmpqKgoKDStgCQlZUFS0tLxdfz5s1Dbm4u9u/fX6NRwnLlRX1hYaHSsfK1Vu9ycxR6e1hMERERUZ2g6sOpXC4HAKxevVrxIfyPanuLak1NzWrbzJo1C5GRkaJYeno6rK2tYWFhgbZt20ImkyE/Px8XLlzAhg0boKGhAXd3dyQmJmLatGm4du2aYv3U2yKXyyEIAuLi4iq9j+qmskkkEvTo0QNJSUmQy+WQyWQYP/7VVvlSqVQxmnf69GlIpVIIgqC4NvBqCuX8+fNrlGu7du3w9ddfV9rmj4WWt7c3Dh8+jNWrV8PT0xMmJibVXgf438/UgwcPRBt+lMcA1VMA6c+HxRQRERHVWa1btwYA6Ovro3fv3lW2tba2BgDcvn1b6ZiqmCotW7ZEfHw8srOz0aBBA9GxmzdvwsjICKampjXqq6J58+Zh7NixolijRo0U/y6VSrF582YcP34cRUVF6NWrFwCgV69eCAwMRFxcHORyebXrpczMzGBgYIBbt24pHbt586ZSrHXr1oiPj4eVlZVS0fA6pFIpYmNjER0djQcPHojyX79+PWJiYpCdnS3K38zMDA0aNMDz58+r/d6W5/rkyRNIpdIqN4/4Y16TJk3CgAED0LNnTyQkJIhGFSvTuXNnAMC5c+eUcjt//jyMjIxUjqTSnw/XTBEREVGd5enpCXNzc4SFheH3339XOl5QUIDc3FwAQKdOndC4cWNERETg2bNnijZ5eXnYunVrja7n7e2NsrIyhIWFieJxcXG4evUqBg0aVOMP8hW1bdsWvXv3Fr0qPuNKKpWirKwMS5YsgZWVFWxsbBTxwsJCrFy5EhKJBG5ublVeR1NTEwMGDMDly5dFO+zJ5XKV27WPGzcOALBw4UKUlpYqHa9uil/F/AEgODgY2tra6NatG4BXD0bW1NRUbP1esZjS0NDAmDFjcPHixUq3Ya+4nfv48ePx22+/VToyVVmuHh4eiI+PR0ZGBnr27Inffvut2vtxd3dH48aNsWPHDuTl5Sni//3vf5GcnIxhw4b9pR/2XJdwZIqIiIjqLH19fezZswfe3t6wtbXFpEmT0KpVK2RnZ+PWrVuIiYnBkSNH4OHhAYlEgjVr1mDMmDFwdnbG5MmTIZFIEBERARMTE6SnpyummFXG19cXkZGR+Oqrr5CRkQE3Nzf89NNP2Lx5MywsLEQ73r1NHh4e0NDQQGpqKnx9fRXxtm3bolGjRrh58yZcXV2r3FWv3PLlyxEXF4cBAwbA398fzZo1w/Hjx/HkyROltp07d0ZISAhCQkLg6OiIYcOGoUmTJnj06BFSUlIQGxuLoqKiaq/p6OiIhg0bIjU1FR4eHopC0cjICJ06dcKFCxfQuHFjpdGv0NBQnDlzBsOHD8fw4cPh6uqKevXq4d69e4iNjYWTk5NiN79Zs2bhu+++w9y5cyGTySCVSmFkZIT79+8jMTEROjo6Slu0l+vevTtOnToFLy8veHh4QCaTVTk9VEtLC9988w1GjBiBHj16wM/PD8+fP8e6detgZmam2Mae/vxYTBEREVGd5unpiUuXLiEsLAz79u3DkydPYGxsDBsbGwQEBIi2PR89ejS0tLSwbNkyBAcHw8LCApMnT0b79u3h4+MDXV3dKq+lpaWFkydPYvny5YiKikJMTAwaNGiAYcOGYfny5aLNDd4mY2NjODo64sqVK0pT+aRSKQ4cOFDtFL9yNjY2+P777zFnzhxs2LAB2tra6NevH/bu3QsLCwul9sHBwejUqRPCw8Oxfv16vHjxAubm5vjwww8RHh5eo2sKggAPDw/ExMQo5dmrVy9cuHBB5S6E9evXx5kzZ7B27VocOnQIx44dg0QiQbNmzdC9e3dMmTJF0VZLSwsnTpzA5s2bsXfvXsVoV5MmTeDs7Kz04OA/cnV1RUJCAvr27Qt3d3fIZLIqv5/Dhg2Drq4uli9fjsDAQGhra6NXr1746quvuF7qL0QoX5ynSlrEnsoP1kDjsJUAAN1bqW/STd1Q/nT0Cn8tIiIiepdSU1PfaB0L/c/atWsRGBiIc+fOwdXVtbbTIaIaquHvwUqHnLlmioiIiKiGioqKlNb+5OXlYdOmTTAxMUHHjh1rKTMiqg2c5kdERERUQ2lpaejXrx9GjhyJFi1a4NGjR4iMjER6ejq2bNmCevXq1XaKRPQesZgiIiIiqiEzMzO4urpi//79ePz4MSQSCdq1a4ewsDAMHz68ttMjoveMxRQRERFRDZmYmODgwYO1nQYR/UlwzRQREREREZEaWEwRERERERGpgcUUERERERGRGlhMERERERERqYHFFBERERERkRpYTBEREREREamBxRQREREREZEaWEwRERERERGpgcUUERER0TuQnJwMQRAQERFR26moLSQkBIIgICMjo0btBUGAr6/vO83pffPw8IC1tbUo5uvrC0EQaiehGvo7fi/+jFhMEREREdFfQmlpKZo2bQpBELBs2bLaTqdGLl68iM8++wzdunWDgYFBlQV2RkYGBEFQ+frwww+V2oeEhODo0aPv+A7Enj17hm+++QZ9+/aFpaUldHV1YWtri08++QS//PKLynNycnLg7++Ppk2bQkdHB/b29tiyZQvkcrmo3Z07d7B48WK4urrCzMwMhoaGcHR0RGhoKF68eKGy79u3b8Pb2xvGxsbQ19dHjx49IJPJ3vp9V0by3q5EREREfxllKcm1nYKIhpNHbafw2tzc3FBQUAAtLa3aTkVtixYtwoIFC6CtrV3bqQAA4uLi8PDhQ9jY2CAiIgKLFi165yNEp06dUvrQ/zpiY2OxadMm2NnZwcHBAWfPnq32nCFDhsDHx0cUa9CggVK7JUuWYMKECfD29lY7v9d14cIFzJkzB7169cLMmTNhamqK69ev4x//+AcOHTqEs2fPom3btor2RUVF6NOnD65evQp/f3+0adMGcXFx+PTTT5GZmYmQkBBF2127dmHTpk0YNGgQxowZAy0tLSQlJWHRokU4dOgQzp8/D11dXUX7n3/+GV27doVEIsG8efNQv359bN++HZ6enoiLi0Pv3r3f+fvBYoqIiIjqnNLSUhQWFkJPT++dXUNDQwM6OjrvrP/3QSKRQCL583xc3LlzJ2xsbPD1119j8ODBSE5ORs+ePd/6dSr+fNSrV++N+po+fTrmzp0LfX19REdH16iYat++PcaOHftG131X7OzscPv2bdjY2IjiH330Efr06YPFixcjOjpaEd+xYwcuXbqE8PBw+Pv7AwD8/Pzw8ccfY8WKFZg4cSKaN28OABg6dCiCgoJQv359xfnTpk1D69atERoaip07d2LmzJmKY0FBQcjOzkZKSgocHR0BAOPHj4e9vT1mzJiBW7duvfNim9P8iIiI6G8tIiICgiAgISEBy5Ytg42NDXR0dHDo0CFFm8uXL2PIkCEwNTWFtrY2bG1tERoaipKSEqX+Dh8+DAcHB+jo6MDKygpLlixBQkKC0vStytZMvXjxAkFBQbCxsYG2tjYaNWqE8ePH4969e6J2Fc/fvXs37O3toa2tjebNm2PVqlU1uvd79+5BEAQEBweL4p6enhAEAevWrRPFXVxc0KZNG8XXla2ZunHjBry8vKCvr4+GDRtizJgxePz4caV5REVFoXv37jA0NISenh5cXFxEH7hrIjMzE99++y3Gjx+P/v37w9zcHDt37lTZ1traGh4eHrhy5QqkUikMDAzQsGFDTJgwQSnP6n4+VK2Zeh0WFhbQ19d/7fNevnyJ/Px8lcfKpwMCQGRkpGg64B+dO3cO7u7u0NfXh4mJCaZMmYK8vLzXzqectbW1UiEFAL1790bDhg1x/fp1UfzAgQPQ09ODn5+fKP7555+juLgYUVFRilinTp1EhVS5ESNGAICo7xcvXuDf//43PDw8FIUUABgYGGDKlCm4c+cOLl26pN5NvgYWU0RERFQnBAYG4p///Cf8/PzwzTffwNbWFgBw4sQJdOvWDXfu3MGcOXMQHh6OLl26YPHixRg1apSoj6ioKAwbNgz5+fkIDg6Gv78//vWvf2HBggU1yqG4uBienp4ICwtDx44dsW7dOowaNQqHDh2Ci4sLfv31V6Vztm7diqVLl2LUqFFYu3YtGjdujPnz5+PAgQPVXq958+Zo2bKlaA1JUVERfvjhB2hoaIjiz58/R0pKCqRSaZV9pqeno0ePHvj+++8xc+ZMLF26FFlZWfDy8lLZftGiRRg5ciQMDQ2xbNkyhIWFQU9PD8OGDcOmTZuqvYdye/bsQWlpKcaPHw+JRIIxY8YgJiYGOTk5Ktv/+uuv6NWrF1q2bIlVq1bBx8cHe/fuRc+ePVUWKZX9fNSGtWvXQk9PD/r6+rC0tMTixYtRWFioOG5mZoa9e/cCAHr06IG9e/cqXhX9+OOPGDBgADp37oyvv/4affv2xc6dOxEQEPDWc87JyUFubi4sLCwUsbKyMly5cgUdOnRQGqV1dnaGIAg1KnjK/7uo2Pe1a9dQWFiILl26KLV3dXUFgPdSTP15xm2JiIiI3qGCggJcvXpVNLXv5cuXmDx5MlxcXCCTyRRT2qZOnQoHBwcEBAQgOTkZHh4eKCkpQUBAAMzMzHDx4kUYGxsDeDWNq3379jXKISIiAmfOnMHcuXNFo0u9e/fGgAEDEBQUpPSB+P79+0hNTVX8xX7SpElo3rw5NmzYgNGjR1d7TalUisjISOTn50NPTw/nz59Hfn4+xo4di2PHjqGkpAQSiQSnT59GaWlptcXUF198gWfPnkEmkymm2M2YMQM+Pj64evWqqO2VK1cQGhqKoKAgrFixQhH/7LPP4O3tjaCgIIwfPx6GhobV3seuXbvg5uamGCWaMGEC1q1bhwMHDmD69OlK7X/++WesW7cOn3/+uSJmb2+PgIAAhIeHKxXAqn4+3jcNDQ1IpVJ4e3ujefPmePLkCQ4dOoRly5bh3LlziI+Ph6amJvT19TF27FiMGzcOLVu2rHRK4LVr13Du3Dm4uLgAePVz/fz5c+zevRtff/01DAwM3lruoaGhKC4uxoQJExSxZ8+eoaCgAE2bNlVqr62tDVNTUzx48KDKfktLS7Fs2TJIJBLRz/vDhw8BQGXf5bHq+n4bODJFREREdcL06dOVPih/9913yMzMxMSJE5GdnY2srCzFq3///gBebUAAACkpKXj48CF8fX0VhRTwalrRtGnTapTDkSNHoKGhgaCgIFH8o48+gqOjI44dO4aysjLRsYkTJ4qmPunp6cHV1RV3796t0TWlUimKi4vx/fffAwBkMhnMzc0xa9Ys5ObmKv56n5SUBEEQqlyDVFZWhuPHj6NTp06idoIgYN68eUrt9+/fD0EQMGHCBNF7m5WVhUGDBiE3Nxfnzp2r9h7Onj2LW7duiT6oOzg4wNHREbt27VJ5jpGRET799FNR7NNPP4WRkRGOHDmi1F7Vz8f7ZmVlhcTERPj7+2PQoEGYPHkyTp48CT8/PyQkJOCf//zna/XXpUsXRSFVTiqVoqSkpMbb3ddEdHQ01qxZAy8vL0ycOFERLx8BrGwDEx0dnUqnMpb7/PPPce7cOSxdulQ0WlhV3+WjYNX1/TawmCIiIqI64YMPPlCKpaamAng12mNmZiZ62dnZAXi1Vgd4Nb0NgMrpXzWdEpaeno4mTZqIirFy9vb2yM3NRVZWlijesmVLpbYmJiZ4+vSp4uucnBz89ttvoldpaSkAKEaayqf0lY8odezYEcbGxqK4g4MDGjZsWGn+jx8/Rl5enuK9qajiDm7lUlNTIZfLYWdnp/T+Tp48GcD/3t+q7Ny5E1paWujQoQN++uknxcvT0xOXL1/GtWvXlM5p2bKl0uYR2traaNmyJdLS0pTaq/r5+LP44osvALyakvo6KvvZASD6+XkTsbGxGDNmDJycnBAVFSVat1VenFacoljRy5cvqyxgv/zyS2zcuBGffPKJ0h8gqur75cuXojbvEqf5ERERUZ2g6oNV+ZbXq1evFi1ir6hJkybvNK/qaGpqVttm1qxZiIyMFMXS09NhbW0NCwsLtG3bFjKZDPn5+bhw4QI2bNgADQ0NuLu7IzExEdOmTcO1a9cwe/bst5q7XC6HIAiIi4ur9D7s7e2r7CMvLw+HDh1CcXExOnTooLLNrl27sH79+jfKtbZHpapiaWkJTU1NpUK7OlX97LzJdu/l4uPj4ePjA3t7e5w6dQpGRkai48bGxtDV1VU53a6wsBBZWVlwd3dX2XdISAiWL1+OiRMnYuvWrUrHy/+7VNV3eUzVFMC3jcUUERER1VmtW7cGAOjr61f7TJrytTq3b99WOqYqpkrLli0RHx+P7OxspecG3bx5E0ZGRjA1Na1RXxXNmzdPad1Mo0aNFP8ulUqxefNmHD9+HEVFRejVqxcAoFevXggMDERcXBzkcnm166XMzMxgYGCAW7duKR27efOmUqx169aIj4+HlZWVaJfA13Ho0CHk5eVhxYoViu9XReHh4di3bx9WrVolGolKS0tDUVGRKFZYWIi0tDSVI2t/ZmlpaSgtLRVtwFDb4uPj4e3tDTs7OyQkJKgcbdXQ0EDHjh1x9epVFBYWiqbkXbx4EXK5HJ06dVI6LyQkRPEMrR07dqjcpbBdu3bQ1tZWOU30/PnzAKCy77eN0/yIiIiozvL09IS5uTnCwsLw+++/Kx0vKChAbm4ugFcfzBo3boyIiAg8e/ZM0SYvL0/lX85V8fb2RllZGcLCwkTxuLg4XL16FYMGDYKGxut/PGvbti169+4telXcPU0qlaKsrAxLliyBlZWVYmtrqVSKwsJCrFy5EhKJBG5ublVeR1NTEwMGDMDly5eRlJSkiMvlcpXbtY8bNw4AsHDhQsW0w4pqOsWvYcOGmDt3LoYOHar0mjx5Mp4+fYpjx46Jznv+/Dk2b94sim3evBnPnz9/rw+5fR2qpt6VlZVh0aJFAICBAweKjhkYGKj8uX3XTp06hSFDhsDW1haJiYlVTg0dNWoU8vPzsW3bNlF8/fr1kEgkim3Pyy1duhRLlizBuHHjsGvXrkr/ezAwMMDAgQORnJyM//73v4p4Xl4eduzYgdatW8PZ2fkN7rJmODJFREREdZa+vj727NkDb29v2NraYtKkSWjVqhWys7Nx69YtxMTE4MiRI/Dw8IBEIsGaNWswZswYODs7Y/LkyZBIJIiIiICJiQnS09OrfUCor68vIiMj8dVXXyEjIwNubm746aefsHnzZlhYWIh2vHubPDw8oKGhgdTUVPj6+iribdu2RaNGjXDz5k24urrWaFe95cuXIy4uDgMGDIC/vz+aNWuG48eP48mTJ0ptO3fujJCQEISEhMDR0RHDhg1DkyZN8OjRI6SkpCA2NhZFRUWVXuvWrVs4e/YsfH19K3148KBBg6ClpYWdO3di2LBhiriNjQ2WLFmC69evw8nJCSkpKdi1axfs7Ozw2WefVXufb8u9e/cUOzTeuHEDAHD8+HHFdt/jxo1TPLTWz88Pz58/R9euXWFpaYmsrCwcPnwYKSkpGDx4MIYOHSrq29XVFQkJCfjqq69gZWUFQRAwcuTI184xIyMDLVq0gLu7O5KTk6tse/nyZQwePBhyuRwTJ05EXFycUpuKo6R+fn7YvXs3AgICkJGRgTZt2iA2NhZHjhzBokWLRM/w2rRpE4KDg2FlZYXevXsrbf9vYWGBPn36KL5euXIlEhMT0bdvX8yePRtGRkbYvn07Hjx4gBMnTrzzB/YCLKaIiIiojvP09MSlS5cQFhaGffv24cmTJzA2NoaNjQ0CAgJE256PHj0aWlpaWLZsGYKDg2FhYYHJkyejffv28PHxga6ubpXX0tLSwsmTJ7F8+XJERUUhJiYGDRo0wLBhw7B8+XJYWlq+k3s0NjaGo6Oj4iG2FUmlUhw4cKDaKX7lbGxs8P3332POnDnYsGEDtLW10a9fP+zdu1flNLTg4GB06tQJ4eHhWL9+PV68eAFzc3N8+OGHCA8Pr/Ja5Q/l9fHxqfLeevbsie+++w6//PKL4j1s1qwZDh06hMDAQBw8eBD16tXDmDFjsGbNGrUeoquu9PR0fPnll6JYTEwMYmJiAADdu3dXFFMfffQR9u7di23btuH333+HtrY27O3tsWnTJkybNk1plGbz5s2YMWMGQkNDFSOo6hRT5efWZI3R9evXFRs8VLbGrmIxVa9ePSQkJGDRokU4ePAgnj59ChsbG2zYsAEzZswQnVe+s+T9+/dFOzeWc3d3FxVTrVq1wpkzZ7BgwQKEhYWhqKgIHTt2RHx8fLXTdt8WoarFZ2kRe95oZVrjsJUAAN1bqW/STd1Q/nT0Cn8tIiIiepdSU1PVXsdCYmvXrkVgYCDOnTuneGAo1R5ra2tYW1tXO8pCr4SHhyMwMBDXr1//U+9q+C7U8PdgpUNcXDNFREREVENFRUVKa3/y8vKwadMmmJiYoGPHjrWUGZH6Tp48ialTp9a5Qupt4DQ/IiIiohpKS0tDv379MHLkSLRo0QKPHj1CZGQk0tPTsWXLFqXnGhH9Fbzu86vof1hMEREREdWQmZkZXF1dsX//fjx+/BgSiQTt2rVDWFgYhg8fXtvpEdF7xmKKiIiIqIZMTExw8ODB2k6DqpGRkVHbKVAdwTVTREREREREamAxRUREREREpAYWU0RERERERGpgMUVERFSHVfW8SSKiv7O38fuPxRQREVEdpaWlhYKCgtpOg4ioVhQUFEBLS+uN+mAxRUREVEeZm5vjwYMHyM/P5wgVEdUZcrkc+fn5ePDgAczNzd+oL26NTkREVEcZGRkBAB4+fIji4uJazoaI6P3R0tKChYWF4vegulhMERER1WFGRkZv/GGCiKiu4jQ/IiIiIiIiNbCYIiIiIiIiUgOLKSIiIiIiIjWwmCIiIiIiIlIDiykiIiIiIiI1sJgiIiIiIiJSA4spIiIiIiIiNbCYIiIiIiIiUgOLKSIiIiIiIjWwmCIiIiIiIlIDiykiIiIiIiI1sJgiIiIiIiJSA4spIiIiIiIiNbCYIiIiIiIiUgOLKSIiIiIiIjWwmCIiIiIiIlIDiykiIiIiIiI1sJgiIiIiIiJSA4spIiIiIiIiNbCYIiIiIiIiUgOLKSIiIiIiIjWwmCIiIiIiIlIDiykiIiIiIiI1sJgiIiIiIiJSA4spIiIiIiIiNbCYIiIiIiIiUgOLKSIiIiIiIjWwmCIiIiIiIlIDiykiIiIiIiI1sJgiIiIiIiJSA4spIiIiIiIiNbCYIiIiIiIiUoPkfVwkLWLPW+mnpe/4t9IPERERERHRm+LIFBERERERkRpYTBEREREREamBxRQREREREZEaWEwRERERERGpgcUUERERERGRGlhMERERERERqYHFFBERERERkRpYTBEREREREamBxRQREREREZEaWEwRERERERGpgcUUERERERGRGlhMERERERERqYHFFBERERERkRpYTBEREREREamBxRQREREREZEaWEwRERERERGpgcUUERERERGRGlhMERER1LXUOQAAD89JREFUERERqYHFFBERERERkRpYTBEREREREamBxRQREREREZEaWEwRERERERGpgcUUERER/X979x98WV3Xcfz5ZncRFsZfsxgJIgo4QBKSWoAFK4qElKMzpZlACEhF0WA6Y1FNYj+txFJUwBQM0YwAzVosGNvhh24gCGPBCgssMBAKIuACuyD77o/P+c5++fb9rncP557PPXefj5k73+XwOfd8zmvuuee+7+ecz5UktWAxJUmSJEktWExJkiRJUgsWU5IkSZLUgsWUJEmSJLVgMSVJkiRJLVhMSZIkSVILFlOSJEmS1ILFlCRJkiS1YDElSZIkSS1YTEmSJElSCxZTkiRJktSCxZQkSZIktWAxJUmSJEktWExJkiRJUguLa3dgS9x+3j+M1O6lxx075p5IkiRJ2to5MiVJkiRJLVhMSZIkSVILFlOSJEmS1ILFlCRJkiS1YDElSZIkSS1YTEmSJElSCxZTkiRJktSCxZQkSZIktWAxJUmSJEktLK7dgS7tuO5eADZet7JuR9pYu5ptdt+7di8kSZIkjciRKUmSJElqwWJKkiRJklqwmJIkSZKkFjZ7z9TMPUhtLXpqQyfPI0mSJEmTxpEpSZIkSWrBYkqSJEmSWrCYkiRJkqQWpup3pmY88PVVP7LNsoMO7KEnkiRJkqaVI1OSJEmS1ILFlCRJkiS1YDElSZIkSS1YTEmSJElSCxZTkiRJktSCxZQkSZIktWAxJUmSJEktWExJkiRJUgsWU5IkSZLUgsWUJEmSJLVgMSVJkiRJLVhMSZIkSVILi2t3QJtsXLsarltZuxutbfPK5bW7IEmSJPXGkSlJkiRJasGRKXVm44BH1cCRNUmSJG0ZR6YkSZIkqQWLKUmSJElqwcv8pIaXKUqSJGlLWExtxgNfXzVSu2UHHdjJ9tatuY31I2yzq+1JkiRJas/L/CRJkiSpBYspSZIkSWrBYkqSJEmSWthq75ka9X6ovra33ZrbeurJ042SQ5f3aHWVu/eNSZIkqTZHpiRJkiSpha12ZEqaNk7tLkmS1C9HpiRJkiSpBYspSZIkSWrBy/ym1KRNsDH07TnhhSRJkuZyZEqSJEmSWnBkSpKkyoY+gQw4iYykrZPFlKSJMA0fJiUNl+9BdVmMa6gspiRJ0jNmMSIN19CP35rFeGRmtY23FREnZeY5tfsxLcyzW+bZLfPsnpl2yzy7ZZ7dMs9umWe3piHPoU5AcVLtDkwZ8+yWeXbLPLtnpt0yz26ZZ7fMs1vm2a3B5znUYkqSJEmSqrKYkiRJkqQWhlpMDfrayglknt0yz26ZZ/fMtFvm2S3z7JZ5dss8uzX4PAc5AYUkSZIk1TbUkSlJkiRJqspiSpIkSZJa6LWYiohtIuLdEbE6ItZHxN0R8aGI2KHr9SPirRFxbkTcGBFPRkRGxO5d71NNfeUZEUsi4qyIuC4iHoiIDRFxR0R8ISIOGM/e1dHza3Rl87qc7/Gq7veufz2+RpdvJsuZx2vGs5f96fn1uSQiTouIm5tj/nsRcVFE7N39ntXRQZ6/HxEXRsTtzWts7WbaHt68j17bbCsjYnlX+zIJes7zLyPiaxHx3eb1eXdE/Os0Zdpznudt5r3zlzrbqYr6yjMidh/hfPSOTneugp5fnxERvxER34yIxyPioYj4SkQc2NkOPROZ2dsD+DsggYuBdwFnAE8CXwW26XJ9YCXwOLAKWN2st3uf+zsteQI7AN8APgKcCpwAnA7cBTwBHFY7i6FlOus1ej9w9DyP59fOYkh5Aj+2QI4nAE8B3wGW1M5jQHkGsKJpewlwMvDHTY7fB/atncWE5JnA94DLgAeBtZtpe17zfnk98M1m3eW1MxhwniuBs4D3AMcDpwE3Nc9xTO0sBpjneU37+d5Hd6udxZDypHxmmi/Ho4GHgA3Astp5DCXPpu0nmvb/CZwCvA+4rclyefUsegz9J4CNwEVzlp/SBPSrXa4P7AYsbv59JlNWTPWd5wLP8ePNgbOidh5DzJTyYWBt7f2eljwXeI63N23/unYeQ8oTeHOz7Ow5bV8KPAZcXjuP2nnO5DHr3/+9ueMZ2AV4VvPv9zJlxVTfeS6w/o6Ugv+m2nkMLU+aYqr2fk9Lngusf1CzrQtr5zGkPIFXNM95Kc3Eec3y5wL3AbcwQvE2zkefl/m9nfJt59/OWf5Jysn56C7Xz8y7MvOHrXs7+XrNcwHfBdYDzxuh7RBUybQZKn92RMQW93iyTcJr9MTm79+P0HbS9Znna5u/585umJm3A1cCr4uI3Ubu+WR6xq+vJo+RZOY9mblhi3o4LL3mucD66yjfdE/DOalKns3lVM+OiGm7p7766xPPR0+zBXnOnI8+k00V1az/EPAlYC+g6mX8fR4sr6ZUsdfMXpiZ64Ebmv8/zvWnTe95RsSiiFgWETtHxKuBz1G+CVzRdicmTI3X6C7AOuBhYF1EXDxF96RUPeYj4iWUN+GrMvPbW9TzydRnns9q/j42z/PMLPuZkXo9uTyndKtKns056QURsX9EnAnsw3Sck2q9Ph9uHo9HxGURMfTjfEbV4z0idgTeCtxJuaxt6PrMc5TzUdV7p/ospl4IPLDAN3P3AMsiYtsxrj9tauS5D+Uen/+lHEBHAH/RPKZB35neAfwV8E7gl4GPA0cC/xUR+7XZgQlT+5g/nvLN2TR8Cwj95vk/zd/DZjeKiKVsKqJeNHLPJ5PnlG71nmfzAfV+yqV9N1CO+XOA3+lyO5X0ned9wIeB3wTeAvw58Crgyoh4fYfbqaX28f42ypfP52bmxjFupy995rnQ+SiAQ5v/rHo+WtzjtpZSbhSbz/pZbZ4Y0/rTpkaedwCHA9sCe1KGcZ9D+dZgGi6p7DXTzHznnDb/HBH/QrmX6gxK1kNW7ZiPiEXAccAjwIUj9HUI+szzs8AfAh+IiEeBy4FllIlnls1qO2SeU7pVI8/HKe+Ti4EXA++gfGBdCjza4XZq6DXPzPy9OYu+GBGfoxSpn6BcSjVktY/3EykjOef+qIYD0Weel1Imlzk5Iu6lTHixFPhd4OWztlVNnyNTj7FpqG6u7Wa1Gdf606b3PDPz0cy8PDNXZOZHKN8SHA5cNFqXJ17112hmXglcAbw2IrbfXNsBqJnnEcCuwOczc1reF3rLMzO/D7yeMlvSOcDtlNHopcAHm7aPjNTryeU5pVu955mZTzXnpK9k5tmUy3p3A74aEUu63FYF1V+fmXkr8E/AnhHxsnFuqwfV8oyIfSmXoV2WmXeNYxsV9JZnM//BkcDVlPPPrcCNwAHAzJcAVc9HfRZT91KG/eYLfxfKcOHmKthnuv60qZ5nc7PvxcAREbHHiP2eZNUzbawFFjH8m6hr5nlC83daLvGDnvPMzG9l5gGUb6QPBfbKzOVsOoGubrEPk8RzSreq55mZTwEXUL6tPmSc2+pB9Twba5u/yzbXaABq5un56BlqJpV7LWUE+lDg5Zm5P5tGwaqej/ospq5ttvfTsxdGxHaUaQ+/Meb1p82k5DkzevL8EdtPsknJdC/KZZMPjth+UlXJMyJeAPwicGNmTtP7QpU8M3NNZl6RmWuaRUdSvgW8ekt3YMJ4TunWpOQ5LeekSclz5vK+7/S0vXGpkmdz39AxlHv7vjSObVRSJc+mqLoiM2fuo3oj5fLJfx/H9kbVZzH1Bco88afOWf4uyqUjF8wsiIg95pnRbOT1txK95RkRO803TWpE7EyZOGEdm24QHLI+M31Oc1/P00TEUZQpPi9rZsUZslrH/LHAEuBT7bo9saq/h0bEKZRv/T+cmUO/J+WZ5qmn6y3PiHjefDe3R8QOlFGA/zfL2AD1mecOzYfgucsPoJzjb87M29o+/4Sodby/CdgJOD8zn+zoOSdB9ffPiHgTcBQl2zu7fv4t6susKdvHv7GIjwK/DVxCmbp0H8qsO1cDh83McBIRa4EXZ2a0Wb9pewibhvl/gTID1Ycovz5NZv7pWHayR33lGRGnUg6YSyiTUDwBvAz4NcqlaCdm5qfHua996THTN1Mmmfgy5X6UH1K+4TmaMiL1msy8ZZz72oc+j/lZ69wM7A68sLn3Z2r0/B66gvLavIly0nwD5cd8/w14yzR8MOggz2Mol51A+bHKbSnnGYA7M/P8WW1/kvLBCsoXJj8PfJryngrw0cx8uMv961tfeTbvn2dT7tddA/wAeAllBGBX4PTMfP9YdrJHPeb5CspN/l+k3I/yKLA/ZXbEjcAbMvOq8exlf/o83metcynlWN83M2/uep9q6vn981OU2XlvoEw887OUCWeuBw6v/t6Z/f5i8iLgPcC3KbOA3EP5QLnjnHZrmeeXuEddv2n7fsoHgHkffe730PMEXkmZ3etWyijUE8DdlG8mDq6dw0Az3YdyY+9tTaYbmn9/DNildg5Dy3NW+4ObY/yC2vs+9DyBP6L8Kv265nEtcDKwqHYOE5TnShY+z6yc0/a4zbRNYPfaeQwlT2APyv0nN1F+E+lJytTeXwaOqp3DAPPcGTifct/JI02edwGfAfauncPQ8pzV/kXAU8DVtfd96HkCvw5c1xzvjwPfAk4Dtq+dQ2b2OzIlSZIkSdOiz3umJEmSJGlqWExJkiRJUgsWU5IkSZLUgsWUJEmSJLVgMSVJkiRJLVhMSZIkSVILFlOSJEmS1ILFlCRJkiS1YDElSZIkSS1YTEmSxioilkdERsR751m20OPAzbRbFxHXR8S7I2JxvT2TJG3tPAlJkmr6PLBinuVrFmgXwM7AscAZwD7ASePsoCRJC7GYkiTVdH1mfnZL20XEx4HVwIkR8QeZef/YeihJ0gK8zE+SNDiZ+SiwijJStUfl7kiStlKOTEmSaloaEcvmLNuQmT8YYd2ZIurBjvskSdJIHJmSJNV0OnD/nMcn52m3NCKWRcROEbFfRHwMOAC4JjNv6a+7kiRt4siUJKmmc4AL5yy7b552pzeP2S4GfmscnZIkaRQWU5Kkmm7NzMtHaDdTdC0B9gPeB+wKrJ+vcUSsAs4ccXILSZJasZiSJA3B7KLr0oi4CrgKOAv4ldkNI2IbSsF1Q79dlCRtbbxnSpI0OJn5NeB84G0RcfDM8ojYGVgHbA+san7g95BK3ZQkTTmLKUnSUP0J8BTwgZkFmXkfcDzld6l2bB5X1OqgJGm6WUxJkgYpM9cA/wi8LiJ+btb/+ing+jq9kiRtTSymJElD9mfARmaNTmExJUnqSWRm7T5IktSZiHgAeGNmXlO7L5Kk6ebIlCRpajQz+T23dj8kSVsHiylJ0tTIzI3A3wD/0czkt2ftPkmSppeX+UmSJElSC45MSZIkSVILFlOSJEmS1ILFlCRJkiS1YDElSZIkSS1YTEmSJElSCxZTkiRJktSCxZQkSZIktWAxJUmSJEkt/B8rDn3pr6J8GgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "font = {'family' : 'sans-serif',\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 18}\n",
    "\n",
    "matplotlib.rc('font', **font)\n",
    "\n",
    "def simpleaxis(ax):\n",
    "    ax.spines['top'].set_visible(False)\n",
    "    ax.spines['bottom'].set_visible(False)\n",
    "    ax.spines['right'].set_visible(False)\n",
    "    ax.spines['left'].set_visible(False)\n",
    "\n",
    "c_ = sns.color_palette('Reds')\n",
    "f,ax = plt.subplots(1,1,figsize=(15,5))\n",
    "\n",
    "#IFR0\n",
    "sns.distplot(dfreport['median IFR0'],color = c_[5],ax=ax, kde=False, norm_hist=True, label='region-wide week 0')\n",
    "mu = np.median(dfreport['median IFR0'].values)\n",
    "ax.axvline(mu, color='r',label='median week 0')\n",
    "print(mu)\n",
    "print(dfreport['median IFR0'].quantile(q=[0.025,0.5,0.975]))\n",
    "\n",
    "#IFRt\n",
    "sns.distplot(dfreport['median IFRc'],color = c_[2],ax=ax, kde=False, norm_hist=True, label='region-wide April 15th, 2020')\n",
    "mu = np.median(dfreport['median IFRc'].values)\n",
    "ax.axvline(mu, color='r',alpha=0.5,label='median April 15th, 2020')\n",
    "print(mu)\n",
    "print(dfreport['median IFRc'].quantile(q=[0.025,0.5,0.975]))\n",
    "\n",
    "\n",
    "ax.set_xlim((0.001,0.2))\n",
    "plt.legend()\n",
    "\n",
    "ax.get_yaxis().set_visible(False)\n",
    "ax.set_xlabel('')\n",
    "plt.subplots_adjust(hspace=0.8,bottom=0.2)\n",
    "\n",
    "x = np.array(list(np.arange(0.01,0.2,0.02)))\n",
    "ax.set_xticks(x)\n",
    "ax.set_xlabel(r'IFR$_{t}$')\n",
    "simpleaxis(plt.gca())\n",
    "plt.savefig(productpath + 'ifrpastpresent.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
